function  [out,errorMes]  = biasCorrectTheta2_threshold_macro_spec(theta2Hat,Var_u,Cov10_u,Cov01_u,x2,macros,...
    econAct,thresholdLevel,bandwidthT,switch_h0,switch_hx,allowDeSaling,numBootStep2,seedNum,numBootStep1And3)
% This function bias-corrects the estimates of Theta2 using a simple bootstrap procedure.
errorMes = 0;
minObsPerRegime = 15;
nm          = size(macros,1);
nx2         = size(x2,1);
factors     = [macros;x2];

% Resampling
rng(seedNum,'twister')
[nx,T]   = size(factors);
% resampling from 2,3,...,T. That is we exclude the first observation
% because here residuals, cov10, cov01 are missing
%indexBoot         = randi(T-1,numBootStep2,T)+1;

% Unfold theta2Hat
[h0Hat_1,h0Hat_2,hxHat_1,hxHat_2,~,gama0Hat,gamazHat,gamaxHat,~] ...
    = unfoldTheta2_threshold_macro(theta2Hat,nx);
residualsx    = zeros(nx,T);
residualsz    = zeros(1,T);
for t=2:T
    if econAct(t-1,1) >= thresholdLevel
        residualsx(:,t) =  factors(:,t) - (h0Hat_1 + hxHat_1*factors(:,t-1));
    elseif econAct(t-1,1) < thresholdLevel
        residualsx(:,t) =  factors(:,t) - (h0Hat_2 + hxHat_2*factors(:,t-1));
    end
    residualsz(1,t) = econAct(t,1) - (gama0Hat + gamazHat*econAct(t-1,1) + gamaxHat'*factors(:,t-1));
end
if all(all(hxHat_1 == hxHat_2)) && all(h0Hat_1 == h0Hat_2)
    % no regime switching
    minObsPerRegime = 0;
end
theta2HatValues     = struc2values(theta2Hat);
theta2BootValues    = zeros(size(theta2HatValues,1),numBootStep2);
AVarTheta2BootDraws = zeros(size(theta2HatValues,1),size(theta2HatValues,1),numBootStep1And3);

s = 0;
while s < numBootStep2
    % resampling from 2,3,...,T. That is we exclude the first observation
    % because here residuals, cov10, cov01 are missing
    indexBoot         = randi(T-1,1,T)+1;
    
    % Generating the data for bootstap s
    factors_Boot     = zeros(nx,T);
    z_Boot           = zeros(1,T);
    Var_u_Boot       = zeros(nx2,nx2,T);
    Cov10_u_Boot     = zeros(nx2,nx2,T);
    Cov01_u_Boot     = zeros(nx2,nx2,T);
    factors_Boot(:,1)= factors(:,indexBoot(1,1));
    z_Boot(:,1)      = econAct(indexBoot(1,1),1);
    Var_u_Boot(:,:,1)= Var_u(:,:,indexBoot(1,1));
    for t=2:T
        if z_Boot(1,t-1) >= thresholdLevel
            factors_Boot(:,t)   = h0Hat_1 + hxHat_1*factors_Boot(:,t-1) + residualsx(:,indexBoot(1,t));
        elseif z_Boot(1,t-1) < thresholdLevel
            factors_Boot(:,t)   = h0Hat_2 + hxHat_2*factors_Boot(:,t-1) + residualsx(:,indexBoot(1,t));
        end
        z_Boot(1,t) = gama0Hat + gamazHat*z_Boot(1,t-1) + gamaxHat'*factors_Boot(:,t-1) + residualsz(:,indexBoot(1,t));
        Var_u_Boot(:,:,t)   = Var_u(:,:,indexBoot(1,t));
        Cov10_u_Boot(:,:,t) = Cov10_u(:,:,indexBoot(1,t));
        Cov01_u_Boot(:,:,t) = Cov01_u(:,:,indexBoot(1,t));
    end
    % Test if we have both booms and recressions in z_Boot
    regime_1       = z_Boot >= thresholdLevel; %a boom
    regime_2       = z_Boot < thresholdLevel;  %a recression
    if sum(regime_1) >= minObsPerRegime && sum(regime_2) >= minObsPerRegime
        s = s + 1;
        % Using the estimator on the bootstrap sample
        theta2Boot = Var1Threshold_GMM_closedForm_macro_spec(...
            Var_u_Boot,Cov10_u_Boot,Cov01_u_Boot,...
            factors_Boot(nm+1:end,:),factors_Boot(1:nm,:),z_Boot',...
            thresholdLevel,allowDeSaling,switch_h0,switch_hx);        
        theta2BootValues(:,s) = struct2array(theta2Boot)';
        if s <= numBootStep1And3
            % We compute the co-variance matrix of theta2Boot for the bootstrap in step 3 of the SR approach
            AVarTheta2Step2Boot = getAVarTheta2_threshold_macro_spec(theta2Boot,Var_u_Boot,...
                Cov10_u_Boot,Cov01_u_Boot,...
                factors_Boot(nm+1:end,:),factors_Boot(1:nm,:),z_Boot',bandwidthT,thresholdLevel,switch_h0,switch_hx);
            AVarTheta2BootDraws(:,:,s) = AVarTheta2Step2Boot.VarRobust;
        end
    end
end

% The bias adjustment - ensuring stationarity of the adjusted estimate
theta2BootAvgValues = mean(theta2BootValues,2);
theta2BootAvg       = values2struct(theta2BootAvgValues,fieldnames(theta2Hat));
biasTheta2Values    = theta2BootAvgValues-theta2HatValues;
biasTheta2          = values2struct(biasTheta2Values,fieldnames(theta2Hat));
%This corresponds to: theta2BiasAdjValues = 2*theta2HatValues-theta2BootAvgValues;
theta2AdjValues     = theta2HatValues-biasTheta2Values; 
theta2Adj           = values2struct(theta2AdjValues,fieldnames(theta2Hat));

% Scaling to ensure stationarity of the bias-corrected model
[params,theta2AdjScaled] = runInduceStatVAR_threshold_macro_spec(theta2Adj,...
    Var_u,Cov10_u,Cov01_u,factors(nm+1:end,:),factors(1:nm,:),econAct,...
    thresholdLevel,switch_h0,switch_hx);

% Co-variance matrix of theta2 from bootstrap
VarBoot = (repmat(theta2BootAvgValues,1,numBootStep2)-theta2BootValues)...
    *(repmat(theta2BootAvgValues,1,numBootStep2)-theta2BootValues)'/numBootStep2;

% Re-centering the draws to account for the bias-adjustment
theta2BootDraws = (theta2BootValues - repmat(theta2BootAvgValues,1,numBootStep2)) ...
    + repmat(struct2array(theta2Adj)',1,numBootStep2);

% The t-tests (under homoskedasticity)
numTheta2          = size(theta2BootValues,1);
tValuesTheta2Asym  = zeros(numBootStep1And3,numTheta2);
for k=1:numBootStep1And3
    tValuesTheta2Asym(k,:)  = (theta2BootDraws(:,k)' - struct2array(theta2Adj))./sqrt(diag(AVarTheta2BootDraws(:,:,k)))';
end
prctilesValues          = [1.0 2.5 5.0 10 50 90 95 97.5 99];
theta2_prctiles         = zeros(numTheta2,size(prctilesValues,2)); 
tTestTheta2_prctiles    = zeros(numTheta2,size(prctilesValues,2)); 
AbstTestTheta2_prctiles = zeros(numTheta2,size(prctilesValues,2)); 
for i=1:numTheta2
    theta2_prctiles(i,:)         = prctile(theta2BootDraws(i,:)',prctilesValues);
    tTestTheta2_prctiles(i,:)    = prctile(tValuesTheta2Asym(:,i),prctilesValues);
    AbstTestTheta2_prctiles(i,:) = prctile(abs(tValuesTheta2Asym(:,i)),prctilesValues);
end
% Storing the percentiles in a struct
for j=1:size(prctilesValues,2)
    label = num2str(prctilesValues(1,j));
    label = strrep(label,'.','_');
    theta2_pTile.(['pct',label]) = theta2_prctiles(:,j);
    tTestTheta2_pTile.(['pct',label])    = tTestTheta2_prctiles(:,j);
    AbstTestTheta2_pTile.(['pct',label]) = AbstTestTheta2_prctiles(:,j);
end
   
%Saving output
out.theta2Hat               = theta2Hat;
out.theta2BootAvg           = theta2BootAvg;
out.biasTheta2              = biasTheta2;
out.theta2Adj               = theta2AdjScaled;
out.delta                   = params;
out.VarBoot                 = VarBoot;
out.theta2BootDraws         = theta2BootDraws;
out.AVarTheta2BootDraws     = AVarTheta2BootDraws;
out.prctilesValues          = prctilesValues;
out.theta2_pTile            = theta2_pTile;
out.tTestTheta2_pTile       = tTestTheta2_pTile;
out.AbstTestTheta2_pTile    = AbstTestTheta2_pTile;
end

